Ciência de Dados I

O objetivo desse notebook é aplicar a Linguagem R em Ciência de Dados na resolução de problemas de forma não supervisionada.

Agrupando dados

A base de dados USArrests contém estatísticas por 100.000 habitantes de prisões por agressão, assassinato e estupro em cada um dos 50 estados dos Estados Unidos em 1973.

data('USArrests')
usa = USArrests
head(usa)

É muito importante tratar valores ausêntes. A função is.na nos diz se um elemento é NA. Combinada com a função any, podemos descobrir se algum elemento da base é NA.

head(is.na(usa))

##            Murder Assault UrbanPop  Rape
## Alabama     FALSE   FALSE    FALSE FALSE
## Alaska      FALSE   FALSE    FALSE FALSE
## Arizona     FALSE   FALSE    FALSE FALSE
## Arkansas    FALSE   FALSE    FALSE FALSE
## California  FALSE   FALSE    FALSE FALSE
## Colorado    FALSE   FALSE    FALSE FALSE

print(any(is.na(usa)))

## [1] FALSE

Também é muito importante normalizar os dados. Mas poderiamos usar normalização por maxmin.

maxmin <- function(x) {
  return (x - min(x))/(max(x) - min(x))
}

normalize <- function(df) {
  return (sapply(df, maxmin))
}

usa = normalize(usa)
head(usa)

##      Murder Assault UrbanPop Rape
## [1,]   12.4     191       26 13.9
## [2,]    9.2     218       16 37.2
## [3,]    7.3     249       48 23.7
## [4,]    8.0     145       18 12.2
## [5,]    8.2     231       59 33.3
## [6,]    7.1     159       46 31.4

Aqui vamos usar a normalização pela média e desvio padrão.

usa = scale(usa)
head(usa)

##          Murder   Assault   UrbanPop         Rape
## [1,] 1.24256408 0.7828393 -0.5209066 -0.003416473
## [2,] 0.50786248 1.1068225 -1.2117642  2.484202941
## [3,] 0.07163341 1.4788032  0.9989801  1.042878388
## [4,] 0.23234938 0.2308680 -1.0735927 -0.184916602
## [5,] 0.27826823 1.2628144  1.7589234  2.067820292
## [6,] 0.02571456 0.3988593  0.8608085  1.864967207

Agora nossos dados estão prontos. Vamos aplicar a função de agrupamento KMeans para 2 grupos.

cls = kmeans(usa, centers=2)
print(cls)

## K-means clustering with 2 clusters of sizes 20, 30
## 
## Cluster means:
##      Murder    Assault   UrbanPop       Rape
## 1  1.004934  1.0138274  0.1975853  0.8469650
## 2 -0.669956 -0.6758849 -0.1317235 -0.5646433
## 
## Clustering vector:
##  [1] 1 1 1 2 1 1 2 2 1 1 2 2 1 2 2 2 2 1 2 1 2 1 2 1 1 2 2 1 2 2 1 1 1 2 2 2 2 2
## [39] 2 1 2 1 1 2 2 2 2 2 2 2
## 
## Within cluster sum of squares by cluster:
## [1] 46.74796 56.11445
##  (between_SS / total_SS =  47.5 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"

Analisando os grupos:

cls$cluster

##  [1] 1 1 1 2 1 1 2 2 1 1 2 2 1 2 2 2 2 1 2 1 2 1 2 1 1 2 2 1 2 2 1 1 1 2 2 2 2 2
## [39] 2 1 2 1 1 2 2 2 2 2 2 2

Conseguimso descobrir qual o número de grupos ideal? Uma alternativa é utilizar uma medida de validação de agrupamento como a silhueta. Para isso precisamos carregar a biblioteca cluster.

require('cluster')

## Loading required package: cluster

sil = silhouette(cls$cluster, dist(usa))
print(sil)

##       cluster neighbor sil_width
##  [1,]       1        2 0.3354254
##  [2,]       1        2 0.3254396
##  [3,]       1        2 0.3911248
##  [4,]       2        1 0.1345762
##  [5,]       1        2 0.3774168
##  [6,]       1        2 0.2972633
##  [7,]       2        1 0.5306315
##  [8,]       2        1 0.1743774
##  [9,]       1        2 0.5035739
## [10,]       1        2 0.4044931
## [11,]       2        1 0.4194608
## [12,]       2        1 0.5680837
## [13,]       1        2 0.3354701
## [14,]       2        1 0.4428995
## [15,]       2        1 0.5945178
## [16,]       2        1 0.5362530
## [17,]       2        1 0.3496924
## [18,]       1        2 0.4387476
## [19,]       2        1 0.5625723
## [20,]       1        2 0.4890198
## [21,]       2        1 0.3837284
## [22,]       1        2 0.5165960
## [23,]       2        1 0.5986325
## [24,]       1        2 0.2620747
## [25,]       1        2 0.1134541
## [26,]       2        1 0.5240161
## [27,]       2        1 0.5926288
## [28,]       1        2 0.4299099
## [29,]       2        1 0.5870108
## [30,]       2        1 0.1782450
## [31,]       1        2 0.5222275
## [32,]       1        2 0.3865530
## [33,]       1        2 0.2596933
## [34,]       2        1 0.5143664
## [35,]       2        1 0.3618181
## [36,]       2        1 0.4064074
## [37,]       2        1 0.1942476
## [38,]       2        1 0.5253841
## [39,]       2        1 0.3635329
## [40,]       1        2 0.3650274
## [41,]       2        1 0.5325247
## [42,]       1        2 0.3175580
## [43,]       1        2 0.3471312
## [44,]       2        1 0.4078063
## [45,]       2        1 0.4408083
## [46,]       2        1 0.2622062
## [47,]       2        1 0.3313743
## [48,]       2        1 0.4592799
## [49,]       2        1 0.5891363
## [50,]       2        1 0.4400334
## attr(,"Ordered")
## [1] FALSE
## attr(,"call")
## silhouette.default(x = cls$cluster, dist = dist(usa))
## attr(,"class")
## [1] "silhouette"

Podemos mostrar graficamente o valor da silhueta por amostra.

plot(sil)

Ou usar a média para avaliar se o valor é bom.

mean(sil[,3])

## [1] 0.408489

Podemos variar o número de grupos para tentar encontrar o valor ótimo de grupos baseado na medida silhueta.

cls = sapply(2:10, function(k) {
  aux = kmeans(usa, centers=k)
  sil = silhouette(aux$cluster, dist(usa))
  mean(sil[,3])
})

Podemos plotar os valores da silhueta.

names(cls) = 2:10
plot(names(cls), cls, type = "l")

Usando a biblioteca plotly.

require(plotly)

## Loading required package: plotly

## Loading required package: ggplot2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

data = data.frame(cls, k=2:10)
plot_ly(data, x=~k, y=~cls, type='scatter', mode='lines')

Vamos carregar duas bibliotecas para nos ajudar a visualizar os dados.

require('gridExtra')

## Loading required package: gridExtra

require('factoextra')

## Loading required package: factoextra

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

Vamos visualizar os agrupamentos.

plot = lapply(2:5, function(k) {
  cls = kmeans(usa, centers=k)
  fviz_cluster(cls, data=usa)
})
do.call("grid.arrange", plot)

Comprimindo Áudio

Um experimento interessante é a tarefa de comprimir áudio utilizando o algoritmo de agrupamento KMeans. Para isso vamos precisar das bibliotecas audio e seewave.

require('audio')

## Loading required package: audio

require('seewave')

## Loading required package: seewave

## 
## Attaching package: 'seewave'

## The following object is masked from 'package:plotly':
## 
##     export

Vamos baixar uma música da internet com permissão de reprodução pelo pixabay. Essa música tem o formato MP3.

url = 'https://cdn.pixabay.com/download/audio/2021/04/07/audio_0bed53371b.mp3?filename=blues-vibes-100-bpm-michael-kobrin-3780.mp3'
download.file(url, 'audio.mp3')

Na sequencia, precisamos converter para WAV e pegar uma amostra para podermos demostrar a compressão. Para isso vamos usar o programa ffmpeg.

system('ffmpeg -i audio.mp3 audio.wav')
system('ffmpeg -i audio.wav -ac 1 mono.wav')
system('ffmpeg -i mono.wav -ss 00:00 -to 00:10 out.wav')

Podemos carregar a música.

sound = load.wave("out.wav")
class(sound)

## [1] "audioSample"

E plotar…

plot(sound, type = "l")

Usando plotly…

data = data.frame(as.numeric(sound), x=1:length(sound))
plot_ly(data, x=~x, y=~sound, type='scatter', mode='lines')

Verificando quantos valores únicos temos nessa amostra.

length(unique(as.numeric(sound)))

## [1] 44923

Para aplicar o KMeans, precisamos transformar esses dados.

wave = data.frame(amp=as.numeric(sound))
summary(wave)

##       amp            
##  Min.   :-0.9010590  
##  1st Qu.:-0.0736106  
##  Median : 0.0000000  
##  Mean   :-0.0000319  
##  3rd Qu.: 0.0752258  
##  Max.   : 0.9559326

head(wave)

Agora podemos aplicar o KMeans. Vamos tentar agrupar os sinais de amplitude da música de 10 até 100 grupos.

require('parallel')

## Loading required package: parallel

cls = mclapply(seq(10, 100, by=10), mc.cores=6, function(k) {
  kmeans(wave, centers=k)
})

O resultado é algo como:

cls[[1]]$centers

##            amp
## 1   0.35133066
## 2  -0.01530472
## 3  -0.32868802
## 4  -0.15544307
## 5   0.23119363
## 6  -0.58907968
## 7   0.70802704
## 8   0.04463795
## 9   0.50747112
## 10  0.13277339

head(cls[[1]]$cluster, 100)

##   [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

cmp = cls[[1]]$centers[cls[[1]]$cluster]
plot(cmp, type = "l")

Agora, podemos salvar alguns desses resultados do KMeans para compressão de música para diferentes valores de k e testar.

cmp = cls[[1]]$centers[cls[[1]]$cluster]
plot(cmp, type = "l")

novoSom = audioSample(cmp)
savewav(novoSom, filename='cmpk10.wav')

cmp = cls[[2]]$centers[cls[[2]]$cluster]
plot(cmp, type = "l")

novoSom = audioSample(cmp)
savewav(novoSom, filename='cmpk20.wav')

cmp = cls[[5]]$centers[cls[[5]]$cluster]
plot(cmp, type = "l")

novoSom = audioSample(cmp)
savewav(novoSom, filename='cmpk50.wav')

cmp = cls[[10]]$centers[cls[[10]]$cluster]
plot(cmp, type = "l")

novoSom = audioSample(cmp)
savewav(novoSom, filename='cmpk100.wav')

Comprimindo Imagens

Um experimento interessante é a tarefa de comprimir imagens utilizando o algoritmo de agrupamento KMeans. Para a compressão, o número de cores precisa ser reduzido. A redução é feita através do agrupamento das cores e seleção daquelas mais frequentes. Para isso precisamos carregar a blbioteca imager:

require('imager')

## Loading required package: imager

## Loading required package: magrittr

## 
## Attaching package: 'imager'

## The following object is masked from 'package:magrittr':
## 
##     add

## The following object is masked from 'package:audio':
## 
##     play

## The following object is masked from 'package:plotly':
## 
##     highlight

## The following objects are masked from 'package:stats':
## 
##     convolve, spectrum

## The following object is masked from 'package:graphics':
## 
##     frame

## The following object is masked from 'package:base':
## 
##     save.image

Algumas imagens são pré-carregadas pelo pacote como a imagem de barcos:

plot(boats, axes=FALSE)

Vamos salvar a imagem original em um arquivo e imprimir seu tamanho em disco:

save(boats, file='boats1.zip')
cat('Tamanho da imagem original:', file.info("boats1.zip")$size, 'bytes')

## Tamanho da imagem original: 1499621 bytes

Podemos extrair informações sobre o objeto da seguinte forma:

boats

## Image. Width: 256 pix Height: 384 pix Depth: 1 Colour channels: 3

dim(boats)

## [1] 256 384   1   3

Gerando o data.frame da imagem boats:

image = as.data.frame(boats)
# cabecalho do data.frame
head(image)

Podemos deixar a imagem borrada:

plot(isoblur(boats, 10), axes=FALSE)

Podemos deixar a imagem em escala de cinza:

plot(grayscale(boats), axes=FALSE)

Podemos remover a cor vermelha da imagem:

test_image = image
test_image[test_image$cc == 1, 'value'] = 0
test_image = as.cimg(test_image)

## Warning in as.cimg.data.frame(test_image): Guessing image dimensions from
## maximum coordinate values

plot(test_image, axes=FALSE)

Podemos separar as cores RGB da imagem da seguinte forma:

red = image[image$cc == 1, 'value']
blue = image[image$cc == 2, 'value']
green = image[image$cc == 3, 'value']

Reconstruindo a base de dados em um data.frame:

data = data.frame(red, blue, green)

Aplicando o KMeans com 10 centroides e 500 iterações máximas:

centroides = kmeans(data, centers=10)
summary(centroides)

##              Length Class  Mode   
## cluster      98304  -none- numeric
## centers         30  -none- numeric
## totss            1  -none- numeric
## withinss        10  -none- numeric
## tot.withinss     1  -none- numeric
## betweenss        1  -none- numeric
## size            10  -none- numeric
## iter             1  -none- numeric
## ifault           1  -none- numeric

Adicionando os centroides encontrados em data:

data[,"cluster"] = centroides$cluster
head(data)

O próximo passo é reconstruir a imagem utilizando as cores identificadas como centroides. A imagem será reconstrída com as 10 cores principais:

new_data = do.call('rbind',
  lapply(1:nrow(data), function(i) {
    centroides$centers[data[i,4],]
  })
)

Plotando a nova imagem:

new_image = image
new_image$value = c(new_data[,1], new_data[,2], new_data[,3])
new_image = as.cimg(new_image)

## Warning in as.cimg.data.frame(new_image): Guessing image dimensions from maximum
## coordinate values

plot(new_image, axes=FALSE)

Vamos salvar a nova imagem em um arquivo e imprimir o seu tamanho comprimido em disco:

save(new_image, file='boats2.zip')
cat('Tamanho da imagem comprimida:', file.info("boats2.zip")$size, 'bytes')

## Tamanho da imagem comprimida: 57574 bytes

Antes de continuar com o experimento, vamos construir a função imageGenerator que recebe uma imagem e o número de centroides e retorna uma nova imagem comprimida:

imageGenerator <- function(image, centers) {
  
  image = as.data.frame(image)
  red = image[image$cc == 1, 'value']
  blue = image[image$cc == 2, 'value']
  green = image[image$cc == 3, 'value']
  data = data.frame(red, blue, green)
  
  centroides = kmeans(data, centers)
  data[,"cluster"] = centroides$cluster

  new_data = do.call('rbind',
    lapply(1:nrow(data), function(i) {
      centroides$centers[data[i,4],]
    })
  )
  
  new_image = image
  new_image$value = c(new_data[,1], new_data[,2], new_data[,3])
  as.cimg(new_image)
}

Usando a função imageGenerator podemos gerar diferentes imagens com diferentes números de cores. Segue o resultado do experimento para 10, 32, 64 e 128 cores:

par(mfrow=c(2,2))
plot = sapply(c(10, 32, 64, 128), function(i) {
  plot(imageGenerator(boats, i), axes=FALSE)
})

## Warning in as.cimg.data.frame(new_image): Guessing image dimensions from maximum
## coordinate values

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 4915200)

## Warning in as.cimg.data.frame(new_image): Guessing image dimensions from maximum
## coordinate values

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 4915200)

## Warning in as.cimg.data.frame(new_image): Guessing image dimensions from maximum
## coordinate values

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 4915200)

## Warning in as.cimg.data.frame(new_image): Guessing image dimensions from maximum
## coordinate values